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Abstract 

Pseudopotentials, tight-binding models, and k-p theory have stood for many years as the standard techniques for computing 
electronic states in crystalline solids. Here we present the first new method in decades, which we call atomistic k ■ p theory. In 
its usual formulation, k ■ p theory has the advantage of depending on parameters that are directly related to experimentally 
measured quantities, however it is insensitive to the locations of individual atoms. We construct an atomistic k ■ p theory by 
defining envelope functions on a grid matching the crystal lattice. The model parameters are matrix elements which are obtained 
from experimental results or ab initio wave functions in a simple way. This is in contrast to the other atomistic approaches in 
which parameters are fit to reproduce a desired dispersion and are not expressible in terms of fundamental quantities. This 
fitting is often very difficult. We illustrate our method by constructing a four-band atomistic model for a diamond/zincblende 
crystal and show that it is equivalent to the sp 3 tight-binding model. We can thus directly derive the parameters in the sp 3 
tight-binding model from experimental data. We then take the atomistic limit of the widely used eight-band Kane model and 
compute the band structures for all III-V semiconductors not containing nitrogen or boron using parameters fit to experimental 
data. Our new approach extends k ■ p theory to problems in which atomistic precision is required, such as impurities, alloys, 
polytypes, and interfaces. It also provides a new approach to multiscale modeling by allowing continuum and atomistic k ■ p 
models to be combined in the same system. 
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I. INTRODUCTION 


An electron in the periodic potential of a semiconductor may be described using pseudopotentials[l, 2], tight-binding 
models[3-6], or k-p theory[7]. All three have been applied to semiconductor nanostructures [5, 8-16] in which the 
translational symmetry is broken by heterostructures, an applied potential, or a finite size. Pseudopotentials and 
tight-binding models are inherently atomistic in that they allow, and even require, specification of the locations of 
atoms. In contrast, k ■ p theory provides a Hamiltonian for the coarse grained crystal using parameters that depend 
on the composition and structure of the material. Alloys are treated in the virtual crystal approximation in which the 
parameters specifying the band structure are empirically fit to the observed electronic properties of a material. While 
this precludes a description with atomic scale precision, typically one does not know the exact position of every atom 
in a system anyway. 

The three methods involve tradeoffs in the approximations made and the physical phenomena that they describe 
most accurately. While k ■ p theory is a continuum model, the momentum matrix elements which parameterize it 
depend on the atomic scale structure of the electronic wave functions. This is advantageous in the computation of 
optical properties since the dipole matrix elements depend on the momentum matrix elements of the Bloch functions 
which also determine the band structure. Tight-binding models use atomistic scale wave functions but involve a 
large number of parameters which are determined using complicated fitting procedures[17] such as genetic algorithms 
[5, 18]. Pseudopotentials also require a large number of form factors and must rely on complex fitting procedures, 
especially when strain is involved[19]. Pseudopotentials are atomistic, but smooth out the core wave function which 
results in smaller momentum matrix elements and thus smaller optical matrix elements. By including enough bands, 
any of the three methods can be made to be accurate throughout the Brillouin zone [20-24]. Here we will focus on 
the dispersion around zone center. 

k ■ p theory in the envelope approximation has been used successfully to describe electronic states in a wide variety 
of inhomogeneous semiconductor systems. The electronic wave function is taken to be a sum of Bloch functions, each 
multiplied by a slowly varying envelope function. The effective Hamiltonian for the envelopes consists of material- 
dependent coefficients multiplying derivatives acting on the envelopes. The electronic states are then determined by 
putting the envelopes on a computational grid and using finite difference approximations for derivatives, giving a 
model that is coarse-grained over a size comparable to the grid spacing. 

An interesting question arises: can the computational grid be made small enough to make the model atomistic? 
Finite differences make the model superficially resemble a tight-binding model with hopping between grid sites, 
suggesting a connection between tight-binding and k-p models. In this paper we will develop an atomistic k -p theory 
by constructing it on a grid in which the sites correspond to atomic positions in the crystal lattice. We will show 
that the atomistic limit of a simple four-band k ■ p theory is equivalent to a tight-binding model, thus allowing the 
determination of k-p parameters and tight-binding parameters in terms of each other. This connection may be used 
to derive atomistic models that identically reproduce the long wavelength physics of k ■ p theory. 

We will begin in Sect. II with a discussion of k ■ p theory in the envelope approximation (henceforth simply k ■ p 
theory) on a three-dimensional grid of points, and using finite differences. We generalize this widely used method to 
an arbitrary grid which need not be cartesian or regular. In Sect. Ill we take the grid to be the crystal lattice itself 
and introduce difference operators on a diamond or zincblende grid. In Sect. IV we discuss how matrix elements 
are computed on such an atomistic grid, and examine the new non-zero momentum matrix elements which appear. 
In Sect. V we take the atomistic limit of a simple four-band model without spin-orbit coupling and show that it is 
identical to a tight-binding model. In Sect. VI we examine the non-Hermiticity that can arise (as in the four-band 
model) and in Sect. VII show how it may be resolved using the finite volume method. In Sect. VIII we take the 
atomistic limit of the more realistic (and widely used) eight-band Kane model with spin-orbit coupling. In Sect. IX 
we describe our fitting method and present our numerical fits of the atomistic parameters for the non-nitride III-V 
semiconductors . We conclude with a discussion of some of the unique features and merits of the atomistic limit. 
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II. ENVELOPE THEORY 


We begin with the basics in order to establish notation. The Hamiltonian for a single electron in a semiconductor 
is 


H = H 0 + H so 

H 0 = Y V 0 (r) Y V e (v) 
ZrriQ 

6 -=ish (,T x w) ■ p 


(la) 

(lb) 

(lc) 


where ho(r) is the crystal potential and V e (r) is a possible externally applied potential. In a translational invariant 
system the single electron states may be found using multi-band k ■ p theory by writing the wave function as 


N 


V>k(r) = V a„ k exp(ik • r)u n (r) 


n—1 


( 2 ) 


where u n { r) are the zone-center Bloch functions and a nk are numerical coefficients. The Hamiltonian is then an N x N 
numerical matrix 


fl 2 A ^ 2 

E, rl 7 j — . k Sm.r>. k- (u m \p\u n ) T (um\H q |w n ) Y (w m |- ~ ,, (cr x W) | , u n ) • k 

4 mic 1 


2mo 


rn 0 


( 3 ) 


where N is the number of zone-center Bloch functions included in the basis. We will initially omit the spin-orbit term, 
and restore it later in section VIII. Note that there is no explicit requirement of small k and Eq. 3 is exact except for 
errors introduced by truncating the basis. For any finite value of N the accuracy decreases with increasing k since 
the solution will not be accurately expressed as a finite sum of zone center Bloch functions. 

If translation symmetry is broken due to a heterojunction or an applied potential then the plane waves of Eq. 2 
are replaced with envelope functions, and the wave function becomes 


N 


V’(r) = ^2f n (r)u n {r), 


( 4 ) 


where the u n (r) are Bloch functions, but otherwise quite arbitrary. Substituting if>(r) into the Schrodinger equation, 
and multiplying both sides by u^(r) we obtain 


h 2 


Q"mn (r)V + P m „(r) ■ V + v mn (r) 


2rno 

a m „(r) = u^(r)u n {r) 

Pmn( r) = —M„(r)pit n (r) 
imo 


fn( r) = Ea mn (jc)f n (v) 


Vmn(r) = U* m (jc) 


2mo 


V(r) 


U n ( r) 


(5a) 

(5b) 

(5c) 

(5d) 


where P mn and v mn are simply functions (i.e. the operators they contain act only on the Bloch functions contained 
within them). We may make one re-arrangement which will be useful when we consider Hermiticity, 


h 2 

2toq 


^mn.(r)V Y Prnniy) * ^ Y V rnn (v') 


f n { r) = Ea mn (r)/„(r) 


Pmn( r ) — 2 (^ mn Pnm ) ■ 


(6a) 

(6b) 


Since u m may be chosen as real, P mn is real, and therefore P m n = ~Pnm ■ The above equations give the effective 
Schrodinger equation for the envelope functions in which the Bloch functions appear as parameters. Note that both 
are exact even with an incomplete set of Bloch functions since the envelopes are arbitrary. For example the trivial 
case of a single Bloch function iti(r) = 1 simply gives back the original Schrodinger equation, in which case the 
solution would consist of an envelope with variation over atomic scales. If a sufficiently large Bloch basis is used, 


3 
















FIG. 1: Wigner Seitz cells in a zincblende crystal, Hi is the cell around the type-I atom at (0,0,0) and fi 2 is the cell 
around the type-II atom at (1/4,1/4,1/4). The hexagonal faces are the planes separating nearest neighbors and the 
small triangular faces are the planes separating second nearest neighbors, which are of the same type. If only the 
nearest neighbor planes were included, the cells would be tetrahedra. Taking next nearest neighbors into 
consideration truncates the corners of the tetrahedra, replacing them with (shorter) triangular pyramids. 


however, the wave function is well approximated with slowly varying envelopes. Eq. 5a will be approximate if the /„ 
are constrained, as they are when defined on a grid that imposes a momentum cutoff. 

To obtain numerical solutions for systems that are not amenable to analytic methods requires reduction to a discrete 
system, such as by using a finite basis set or functions defined on a grid. For nanostructures with irregular geometries, 
putting the envelope functions on a grid and replacing derivatives with finite differences is especially convenient since 
no assumptions about the geometric symmetry are required[ L3, 25, 26]. We denote the grid points with coordinates 
R and use r for continuum coordinates. The continuous space can be broken up into cells Hr, each centered on the 
grid site at R, and the integral over all space can then be written as a sum of integrals over the cells 

[d 3 r = J2[ d 3 r. ( 7 ) 


If a sufficiently large number of Bloch functions are used the envelope functions will be slowly varying and /„(r) 
and its derivatives will be approximately constant over each cell. This seemingly reasonable assumption can result 
in a non-Hermitian Hamiltonian, which will be resolved in Sect. VI. Integrating Eq. 5a over Hr and approximating 
derivatives of f n ( r) by finite differences on the grid we obtain 


h 2 

2mo 


(?x m |'u n )o R A“ -T 


(^m |P |^n) C2 r 

wtiq 


('U j rn\H 0 1 U n )o F 


In R 


-^(^m |^n) Or /nR, 


( 8 ) 


where f „r is the nth envelope function on the site at R, and A is the finite difference approximation to the gradient, 
A/„r ~ d x f n (r) | r _ R , which is a weighted sum of the values of /„ at R and nearby grid sites. We adopt an abbreviated 
notation for the projected matrix element of an operator O, 


/ u* m (r)Ou n (r) d 3 r = (u m \d\u n )n R . (9) 

J Qr 

If Hr contains an integer number of crystal unit cells then (M m |M ra )n R = 5 mn . The solution of Eq. 8 is obtained by 
computing the eigenvalues and eigenvectors of a large sparse matrix, for which there are efficient algorithms [27, 28]. 
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If the Bloch functions are the same throughout the structure then the ('u m |p|u n )n R are constants and Eq. 8 can 
be used directly to determine the electronic states. This will be the case if confinement is provided by an externally 
applied potential or for a nanocrystal in which the vacuum is modeled as a large potential barrier. In a heterostructure, 
however, the matrix elements appearing in Eq. 8 will vary spatially. In the atomistic limit, in which the Rs correspond 
to individual atoms, the matrix elements will vary spatially even in a bulk crystal if it contains different atoms. This 
will cause the Hamiltonian in Eq. 8 to be non-Hermitian, requiring a more careful treatment. We will return to this 
problem and its remedy in section VI. 


III. ATOMISTIC GRID 


The finite difference approximation consists of replacing a derivative at a grid site with a difference operator which 
acts by taking weighted sums of the values on nearby grid sites. For example, on a uniform cartesian grid the derivative 
of a function / at a grid site located at R may be approximated using the symmetric difference 


A*/(r) 


r=R 


/(R + ex) — /(R — ex) 
2e 


9 x f( r) 


+ 0(e) 

r=R 


( 10 ) 


where e is the grid spacing and x is a unit vector. This gives the lowest order approximation to d x , and more accurate 
results may be obtained by including more sites in the sum[29]. The accuracy of the finite difference approximation 
improves as the grid spacing decreases, making it tempting to shrink the grid spacing as much as possible. The 
existence of a physical crystal lattice suggests using the crystal lattice itself as the computational grid. The values of 
the envelope functions will then be defined on the atoms themselves, yielding an atomistic theory. We will develop 
this model for the diamond/zincblende lattice because of its importance in semiconductor physics, but the approach 
can be applied to any crystal structure. 

On a rectangular grid the low order finite difference approximations may be written down intuitively, but on a 
non-rectilinear grid one needs a more systematic approach. The general method for constructing a difference operator 
is to write down a Taylor series expansion of a function at a point R, express the function values on sites in terms of 
that expansion, and solve for the linear combination of the function values on the site and its neighbors that gives the 
desired derivative to lowest order [29-32]. This method is used to obtain high-order difference approximations with 
smaller errors, and it may be used with non-rectilinear grids or even irregular grids. 

Because the zincblende/diamond lattice has a basis with two atoms per unit cell the difference operators will be 
different on the inequivalent sites. We denote the atom at (0, 0, 0) as type 1 (assumed to be the anion in zincblende) 
and the atom at 2i j a -(l, 1,1) as type 2 (cation in zincblende), where ai a tt is the lattice constant. The nearest neighbors 
of the type 1 atoms are at displacements 


di 

d 2 

d 3 

d 4 


dlatt 

~T~ 

&latt 
O'latt 
C^latt 


(1,1,1) 

(-i,-i, i) 
(-i,i,-i) 
(i,-i,-i) 


(lla) 

(lib) 

(11c) 

(lid) 


and the nearest neighbors of the type 2 atoms are at displacements —d n . On the type 1 sites the nearest neighbor 
differences are given by 

ds/(r)| r _ R = Aia;/(R) + O(ciiatt) = [/(R + di) — /(R + d 2 ) — /(R + d 3 ) + /(R -I- d 4 )] /cii a tt + 0(ai a tt) (12a) 

d v f(v) | r=R = A ly /(R) + 0(a latt ) = [/(R + di) - /(R + d 2 ) + /(R + d 3 ) - /(R + d 4 )j /a latt + 0(a latt ) (12b) 

d*/(r)j r=R = A lje /(R) + 0(a la tt) = [/(R + d 4 ) + /(R + d 2 ) - /(R + d 3 ) - /(R + d 4 )]/a,„« + 0{a latt ) (12c) 

V 2 /(r)| r=R = A ?/( R ) + 0(a latt ) = 4- [/( R + di) + /(R + d 2 ) + /(R + d 3 ) + /(R + d 4 ) - 4/(R)] + 0(a la tt) 

a latt 

(12d) 
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and on the type 2 sites the difference operators are 


A 2x /(R) = [ - /(R - di) + /(R - d 2 ) + /(R - d 3 ) - /(R - d 4 )] /a lat t (13a) 

A 2 „/(R) = [ - /(R - di) + /(R - d 2 ) - /(R - d 3 ) + /(R - d 4 )] /aiatt (13b) 

A 2z /(R) = [ - /(R - d 4 ) - /(R - d 2 ) + /(R - d 3 ) + /(R - d 4 )] /oz„« (13c) 

Al/(R) = [/(R - d 4 ) + /(R - d 2 ) + /(R - d 3 ) + /(R - d 4 ) - 4/(R)]. (13d) 

®latt 


Using the site at R and its four nearest neighbors, only the four derivatives V and V 2 can be constructed. Substituting 
the above difference operators in Eq. 8 gives the Hamiltonian for the envelopes, parameterized by the matrix elements 
in Eq.s 8 which would be empirically fit to measurements on bulk materials. 

Perturbative k ■ p models require approximations for second derivatives as well. Since no derivative approximation 
beyond V and V 2 can be constructed with four nearest neighbors, second nearest neighbors must be used. The second 
nearest neighbor differences are the same for type 1 and 2 sites, 


d 2 J( R) = 


A latt . 


4/(R) — /(R + d-q-qo) — /(R + di,_qo) — /(R + d_i j4j o) — /(R + dqqo) 
— /(R + d_ 4 i o,-i) — /(R + dqo,-i) — /(R + d_ 1,0,1) — /( R + dqop) 
+/(R + do, _ i,_i) + /( R +do,q_i) + /(R + do,_i,i) + /( R + do,i,i) 

+^( a Lti) 


cW(x) | X=R = 


*latt 


/(R + di,_i,o) + /(R + d_i,i,o) — /(R + dqqo) — /(R + d_i,_i,o) 


0(a 2 


latt) 


(14a) 

(14b) 


where dqyfc = [ix + j y + kz] is the displacement to the second nearest neighbor from the central site at R. 
Difference formulas for d 2 , 9 2 , d x d x , and d y d z are obtained by cyclic permutation of xpy,z. Note that while the 
approximation to V 2 involves only nearest neighbors, mixed derivatives and d xl S 2 , and <9 2 require next nearest 
neighbors. This does not present any fundamental or technical problems and is consistent with the source of these 
terms, second order perturbation theory. Since the zeroth-order Hamiltonian contains nearest neighbor couplings 
from V 2 , the second order Hamiltonian contains next nearest neighbor couplings. 


IV. MATRIX ELEMENTS 


The atomistic matrix elements are somewhat different from those usually appearing in k ■ p theory due to the 
projection of the Bloch states to atomistic cells. The most obvious choice for Hr would be the Wigner Seitz cell 
around each R, as shown in Fig. 1 for the diamond/zincblende lattice. Many of the same selection rules from k ■ p 
theory apply since they depend on symmetry, which the atomic cells possess. The most notable difference is 
the existence of matrix elements that are zero in the continuum model but are non-zero in the atomistic limit with 
cancellations between the atomistic cells. Over the primitive unit cell, Hi + U 2 , the Bloch functions satisfy 


('Urn |^n) O1+Q2 — (^m|^n)fii 4“ (^m|^n)n 2 — ^mn 


(15) 


but the atomistic matrix elements are not necessarily proportional to S mn since for m ^ n there could be two nonzero 
terms that cancel. If no two bands transform as the same representation of Td, then (u m |rt„)n 12 oc S mn . For example 
in a model with an SXYZ basis, (A|5)n 12 = 0 because under a rotation by tt about the y-axis the crystal is invariant, 
but (A|5)n 12 will change sign. In this paper we will consider models in which all the states transform differently, and 
in particular models derived from an SXYZ basis. In these cases the right hand side of Eq. 8 is a diagonal matrix, 
which we choose to be a multiple of the unit matrix, and there is no need to solve a generalized eigenvalue problem. 
It only necessary to multiply the left hand side by the inverse of this matrix. 

The use of atomistic cells modifies the momentum matrix elements as well. For states projected to a volume H, the 
momentum matrix element is given by 


[ u * m ( r ) -V u n { r) d 3 r = 

Jn * 



tV[i 4 (r)«„(r)] err 


(16) 
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where the second integral on the right hand side will vanish due to periodicity if fl contains an integer number of unit 
cells. If S7 contains some fraction of a unit cell, then uj^(r)u n (r) is not periodic over 17 and there is no reason for the 
second integral on the right hand side to vanish. We write the projected matrix element in the more compact form 

(^m|p|^n)ni = ('^n|p|'^m)Oi d“ Ilf linm (17) 


where Tln in m is given by the second integral on the right hand side of Eq. 16 and is nonzero only if 17 does not contain 
an integer number of unit cells. If we sum over two sub-cells 17i and 172 to make a whole unit cell the correction must 
vanish, and therefore n ainm = -IIo 2rim . 

The atomistic momentum matrix elements between the conduction and valence bands obey the same selection rules 
as in continuum k ■ p theory since they rely only on Td symmetry, however because the projected matrix elements on 
the two atoms can be different we have 

iPo = —{S\p x \X)n 1+ n 2 = — (S\p x \X) ni + — (S\p x \X)n 2 = iP ai + iP a2 (18) 

mo m o mo 

where the subscript a denotes that the matrix elements are projected to single atoms. Throughout this paper we will 
use an a subscript to distinguish atomistic parameters from those of the continuum k ■ p theory. For the diamond 
crystal P ai = P a2 due to inversion symmetry, while for zincblende P ai ^ P a2 . 

In k ■ p models with more than one p-like band, such as the 16 and 14-band models[33, 34], there are also matrix 
elements of the form iQ = — (X v \p y \Z c ) where the v and c subscripts indicate the valence and conduction bands. 
Matrix elements with this general form but within the same band, such as (X v \p y \Z v ), are zero by a simple symmetry 
argument which depends on the periodicity of the unit cell. Due to invariance under a 180° rotation about the 
y-axis, ( X v \p y \Z v ) = (Z v \p y \X v ). For matrix elements over a whole unit cell we also have (X v \p y \Z v ) = (Z v \p y \X v )*, 
and therefore ( X v \p y \Z v ) = {X v \p y \Z v )* . But since the states |X), |Y), | Z) can all be taken as real, (X v \p y \Z v ) 
must be pure imaginary and therefore the matrix element is zero. In the atomistic case, we see from Eq. 17 that 
{X v \p y \Z v )Q i (X v \p v \Z v )q. and the above argument breaks down. While the projected matrix element can be 
nonzero, the matrix element over the whole cell is zero and therefore 


0 — (X v \p y | Z v }q 1 -(- (X v \p y | Z v ^q 2 iQa 'iQa- 

mo mo 


(19) 


We see that the atomistic k -p model has an additional parameter not present in the continuum k-p model from which 
it is derived. In an inversion non-symmetric crystal the atomistic limit will also double the number of Hamiltonian 
matrix elements since there will be different matrix elements on each atom. 

Summarizing, in the atomistic model we have 


iP ai = —(S\p x \X) Ui = —(S\p y \Y) Qi = — (S\p z \Z) Qi 

777-0 777-0 777-0 

iP' ai = ~ — (X\p x \S) ni = ~ — {Y\py\S) ni = ~ — {Z\p z \S) ai 

777-0 777-0 777-0 

iP ai + iP a2 = -iP' ai - iP ' a2 = iPo 


(20a) 

(20b) 

(20c) 


where Pq is the usual continuum k • p parameter. In addition, there are new intraband matrix elements 




P(X\p y \Z) Qi = —(Z\p x \Y)a = P(Y\p z \X) ni 

mo mo mo 

— (Z\p v \X) Qi = —(X\p z \Y)n z = P(Y\p x \Z)n, 


( 21 ) 


which satisfy iQ ai = — iQa 2 ■ 

As will be seen in Sec. V, in an inversion non-symmetric crystal the combination of finite differences and the fact 
that Pi P' results in a non-Hermitian Hamiltonian. One solution is to simply use an inversion symmetric basis, 
which is reasonable since k ■ p theory is often formulated in the symmetric approximation. As will be shown in Sect. 
VIII, inversion symmetry may still be broken by the sub-unit cell structure of the envelope function. Alternatively, 
we may change the volumes of 17i and H 2 by using generalized Voronoi cells[35, 36]. By adjusting f7i and H 2 we can 
make Pi = P' while maintaining inversion non-symmetry. A generalized Voronoi cell may be constructed by rescaling 
the distances that would be used to determine the Wigner Seitz cell. Consider a site at R, with nearest neighbors at 
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Raw.i, and next nearest neighbors at R nnnj- The generalized Voronoi cell around R is the set of points p satisfying 
the two conditions 


o p — R| < 

P ~ Ratjv, i 

(22a) 

p — R < 

P “ R-ATAfAT,i 

(22b) 


where a is a scaling factor that determines the relative sizes of a cell. Note that the second condition does not have a 
scaling factor because the next nearest neighbors are of the same type. Using s / 1 on the type-1 sites, and a —> 1/a 
on the type-2 sites we change the relative volumes of Hi and n 2 while maintaining fli + fl 2 as a unit cell. Shrinking Hi 
will decrease both P ai and P' ai while increasing P a2 and P' a2 , and therefore we may adjust the scaling factor to make 
P ai = P' a2 and P' = P a2 . This method will prove useful for restoring Hermiticity in the inversion non-symmetric 
case in Sec. V. 


V. FOUR-BAND MODEL 


To demonstrate the basic structure of the atomistic limit we first consider the simple four-band k ■ p model without 
spin orbit coupling, using zone-center Bloch states |X), |y), \Z) for the valence band and | S) for the conduction 
band. This model and the tight-binding model with which it will be compared are too simple to be used for realistic 
calculations, but they demonstrate the basic structure of the atomistic limit. Using plane waves for the envelopes, 
the k ■ p Hamiltonian is 


H 0 = 


S 


/ E c + £ c (k ) 
I -iP 0 k x 
I -iP 0 k v 
\ -iP 0 k z 


X Y 

iP 0 k x iP 0 k y 

E v + £ v (fc) 0 

0 E v + £ v (k ) 

0 0 


Z 

iP 0 k z \ 
E v + £ v (k)) 


iP 0 

£ c (k) 

£ v {k) 


h 

mo 


(S\p x \X) 



(23a) 

(23b) 

(23c) 

(23d) 


where F c and F v are remote band contributions to the conduction and valence band respectively, and are included 
to make the k ■ p model agree with a tight-binding model. In the atomistic limit the Hamiltonian also includes terms 
involving the momentum matrix element of Eq. 21 


S 

X 

Y 

Z 


O 

O 

0 

0 1 


0 

0 

iQak z 

iQ a^y 

(24a) 

0 iQak z 

0 

i'Qak’X 

V 0 i'Qaky 

iQakx 

0 ) 


(X\Py\Z) Ql = 

~{X\p y \Z)n 2 

too 

(24b) 


Eq. 24a is only for an anion site, and on the cation site we will have —Hq. In a plane wave basis, it is convenient to 
define the approximate wave vectors obtained from the finite difference operators acting on plane waves, 


IC ln — ie 


-ik r / 


^1 n ^ 


ikr 


(n = x,y,z) 


JC( = -e 


-ik r a 2 _ik-r 


^2n = -ie 1 r A 2 „ e 
K\ = —e _ik r A2 e ik '' 


ik 


r = ic ln 
= k\* 


(25a) 

(25b) 

(25c) 

(25d) 


(n = x,y,z 




where the subscripts 1 and 2 on K, indicates the atom type on which the difference is centered. On a dia- 
mond/zincblende lattice, using Eq. 8 and the finite differences defined in Eq. 12a-13d we obtain 


Si 

Xi 

V 

z i 

s 2 

X 2 

y 2 

Z 2 

^ E c i + V c 

0 

0 

0 

Td 

iP ai Pix 

iPa\ Ply 

iP ai Piz ^ 

0 

E v i + V v 

0 

0 

-iP’ ai Pix 

Tvi 

iQa Plz 

iQa Ply 

0 

0 

E v i + V v 

0 

~iP’ax P'4y 

iQa P'1 z 

Tvl 
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where the subscripts 1, 2 label the atoms within the unit cell and 
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where b is the band index ( v or c for valence or conduction), c indicates the type of atom (1 or 2), and the momentum 
matrix elements are given by Eq.s 20c and 21. 

If the crystal lacks inversion symmetry, then P ai ^ , P' ^ P a2 , and H is not Hermitian. The problem may 

be remedied by starting with a set of inversion symmetric Bloch functions, in which case P ai = P ao = P' a = P' a . 
The Hamiltonian can still be inversion non-symmetric due to the different potentials on the anion and cation, giving 
rise to inversion non-symmetric zone-center envelope functions. An alternative approach is to modify the differencing 
scheme using the generalized Voronoi cells described at the end of Section IV. The cell size can be adjusted so as to 
make P a2 = P' a and P' a2 = P ai , restoring Hermiticity while maintaining P ai ^ P 1 and P a2 ^ P' . Deforming the 
cells will also change the values of the Q ai , but since Q ai = —Q a2 the form of H will not be affected. In a model 
with more bands, modifying D to make H Hermitian would appear to be limited to tuning the P ai s for just one pair 
of bands. One could use different Ds for different bands, or simply set the P ai s on an ad hoc basis. 

We may compare the atomistic four-band k ■ p model of Eq. 26 with the four-band tight-binding Hamiltonian[37 
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where the gs are the standard tight-binding functions, which are related to the P s by 
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Equating the matrix elements in Eqs 26 and 28, the tight-binding and k ■ p parameters are related by 
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We thus find that the atomistic limit of the k ■ p model is equivalent to a tight-binding model. In order to make this 
one-to-one correspondence it was necessary to include spherically symmetric remote band contributions to both the 
conduction and valence bands and to have different momentum matrix elements on the two atoms (at least in the 
inversion non-symmetric case). Our inclusion of only limited remote band contributions to the valence band gives the 
Luttinger model in the spherical approximation with 7 = f (272 + 373 ). 

The atomistic limit always adds at least one new parameter, Q a . For inversion symmetric crystals the Hamiltonian 
matrix elements are the same on both atoms, so there are as many parameters as in the original k ■ p model, plus 
the additional parameter Q a . For inversion non-symmetric crystals the Hamiltonian matrix elements are different on 
each atom, so the number of diagonal matrix elements is doubled, plus the additional Q a . In Sect. VIII we will take a 
hybrid approach in which the Hamiltonian is not inversion symmetric, but the Bloch basis is chosen to be symmetric. 
In that case the matrix elements of FI 0 (c.f. Eq. la) will be different on different atoms, but the momentum matrix 
elements will be the same on each atom since they depend on derivatives of the inversion symmetric Bloch functions. 
Therefore, the number of diagonal parameters will be doubled, the number of momentum matrix elements will remain 
the same, and Q a will be added. 

An important feature of the atomistic limit is that the envelope varies within a unit cell even at zone center, and 
thus modifies the effective Bloch functions. In Eq. 26 we see that A 2 couples the two atoms even at k = 0 via 
the off-diagonal matrix elements T c i, T v 1 , T C 2 , and T V 2 ■ This results in a doubling of the number of bands over 
the continuum model, with the additional bands being shifted by an energy on the order of h 2 /moa 2 att . Since the 
atomistic model includes two grid sites per unit cell the envelope functions include wave vectors outside the first 
Brillouin zone. These states may be interpreted as approximate Bloch functions with different symmetry from the 
zone-center Bloch functions of the theory. For example, when multiplied by an envelope that changes sign from site 
to site the anti-bonding S-like Bloch function of the conduction band becomes similar to the bonding S-like state. Of 
course such a ’’fake” Bloch function is not the true zone-center Bloch function, but an approximation. This simple 
model illustrates the basic features of the atomistic limit, but in order to develop more realistic models we need to 
examine the inversion non-symmetric case more closely and study the relationship to heteroj unctions. 


VI. NON HERMITICITY 

As we saw in Sec. V, straight-forward application of the atomistic limit to an inversion non-symmetric crystal gives 
a non-Hermitian Hamiltonian since the momentum matrix elements are different on different atoms. This problem 
generally arises when a finite difference operator is multiplied by a spatially varying coefficient, and also occurs at a 
heteroj unction in continuum k ■ p theory [13, 25, 26, 38]. In this section and Sec. VII we will examine this issue in a 
general framework that is applicable to both continuum k ■ p models that have been put on a grid and our atomistic 
model. Consider a Hamiltonian containing a Hermitian differential operator T> which is approximated by a difference 
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operator D consisting of a weighted sum of the function values on nearby grid sites, 
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( 31 ) 


where (1 r r/ are coefficients defining the difference operator, most of which are zero except when R' and R are close 
to each other. Since D is Hermitian, (g|D/) = (Og|/) and the corresponding difference operator must satisfy 


E 5r ^ R R ' / R ' — E f ^RR' 9b. 1 /r — E ^R'.R 9r /r' (32) 

R.R' R,R / R,R' 

where the sums on R and R' range over all lattice sites. Therefore the Hermiticity of D requires gIr.r' = c!r, r . If the 
continuum Hamiltonian contains a differential operator multiplied by a spatially varying coefficient c(r) then simply 
multiplying that operator by c(R) gives 
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(33) 


Since Hermiticity requires crgIr r/ = [cr/^r/ .r]*, any spatial variation in the magnitude of c(r) spoils the Hermiticity. 
This problem arises in a one-band model with a spatially varying effective mass as well as in multi-band envelope 
models with spatially varying parameters. This problem will occur in Eq. 8 for a heterostructure, but in the atomistic 
limit it will arise even for a bulk material if the atoms differ from one another. 

A common solution is to symmetrize over the connected sites[I3, 25, 26], 
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The symmeterization is applied to A x , A y , and A z . which are then used to construct other operators. This resolution 
of the problem is not as ad hoc as it may seem since the first derivative is most naturally defined on a link between 
two sites. For example on a one dimensional grid with spacing e the difference between two adjacent sites gives an 
approximation to the derivative on the link connecting them, 


3*/(z)L Xo+e/2 ~ A*/0z)U 


xo+e/2 


= - (/(*0 + e) - f(x q)) . 


(35) 


Therefore, the value of a coefficient multiplying the derivative is naturally defined on the link itself and should be 
interpolated between the two points being differenced, giving 


c(x)d x f(x)\ x=xo+e/2 
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2 e 


[c(x 0 ) + c(x 0 + e)] [f(x o + e) - f(x 0 )}. 


(36) 


For terms containing variable coefficients and second derivatives, one must also be careful about the operator ordering. 
Neither c(x)d x nor d x c(x) can be made self-adjoint, even in the continuum, however a symmetrized operator such as 
d x c( x)d x can [38-40]. Therefore we can write 
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(37) 


VII. FINITE VOLUME METHOD 

The symmetrization procedure described above has an intuitive appeal, however a more formal approach will give 
the same result while providing some additional insight into the problem. The root cause of the non-Hermiticity is 
that in discretizing Eq. 5a we assumed the envelopes and their derivatives were constants over Hr. Instead, we can 
make use of the finite volume method [ 11] in which the divergence theorem is used to convert the volume integral over 
a cell into a surface integral. Discretizing this modified version results in finite differences over two sites that are 
multiplied by a quantity defined on the link between the sites. This means that the coefficient is the same (except for 
a possible sign) when evaluated on either of the sites, thus guaranteeing Hermiticity. 


11 











Let us return to the continuum, but using Eq. 6a rather than Eq. 5a, and consider the V mn { r) ■ V/„(r) term. 
When integrated over a volume around the grid site R, 


(38) 


[ P mn (0 ■ V/„(r) d 3 r = f f n (r)V mn (r) ■ ds - [ /„(r)V • V mn (r) d 3 r 
J j c7F2r j $7r 

where SHr is the bounding surface of Hr. Since f n (r) is slowly varying we can replace it with f n r in the second 
integral on the right to obtain 

[ Pmn(r)'V/„(r)d 3 r« / /n(r)7mn(r) ■ ds - f nR f V mn {r) ■ ds. (39) 

J ^R J 9Qr J 5Qr 

The surface <9Hr is a polyhedron centered on the point R with faces Sd , each of which is normal to the displacement 
from R to the neighboring site at R + d (see Fig. 1). Because / is slowly varying, its value on one of the faces Sd 
is approximately constant, with a value that may be approximated by linearly interpolating between the two sites, 
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and Ad is the forward difference in the d direction defined by Ad./n r = [/ n R+d — /„r.] /1 d |. The discretized k-p term 
now consists of a link connecting sites multiplied by a coefficient defined on the link, which is therefore Hermitian. 
When approximated with naive finite differences, V • a mn (r)V/„(r) also becomes non-Hermitian. Applying the same 
methods as above gives 
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Using the finite volume method, a derivative term becomes a sum of differences between sites multiplied by a 
coefficient that depends on an integral over the surface separating the cells centered on the two sites. Therefore 
the matrix element connecting two sites will be the same whether evaluated at R or R + d, making H Hermitian. 
The quantities given by Eq.s 41b and 40c do not need to be explicitly computed and we may simply use the naive 
differencing formulas with coefficients empirically fit to bulk properties. The finite volume method provides the 
justification for the coefficient being determined by the two sites being connected. 

We may see explicitly how the finite volume method works by considering the Hamiltonian matrix elements involving 
the S and X states on two nearest neighbor sites at Ri and R 2 . Using the differences given by Eq.s 12a and 13a in 
Eq. 8 we would have 
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which is clearly not Hermitian since all four of the non-zero matrix elements are different. Using Eq. 6a will partially 
symmetrize Eq. 42 because V mn = —Vum, but the matrix elements projected to Hi and H 2 will still be different if 
the basis is inversion non-symmetric. Applying the finite volume method then gives 
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where we have omitted the coordinate indices on V . Since V mn is real and V m n = —Pnm, H is Hermitian. 

There are actually two separate symmetrizations: Eq. 6a and the finite volume method. The former anti- 
symmetrizes the momentum matrix element with respect to the band index (’’index symmtrization”) and the latter 
symmetrizes with respect to the coordinate (’’spatial symmetrization”). With an inversion symmetric Bloch basis 
only the index symmetrization is necessary to obtain a Hermitian Hamiltonian, but in the inversion non-symmetric 
case the finite volume method provides spatial symmetrization. The finite volume method results in equal momentum 
matrix elements, unlike the generalized Voronoi cell approach used in Sect. V which gave two distinct S-X momentum 
matrix elements for an inversion non-symmetric basis. It is interesting that Eq. 6a gives a diagonal kinetic term with 
the same form as the common operator ordering choice in the effective mass approximation, V • m ( r ) V?/>(r)[38— 40]. 

Note that if the Bloch functions centered on R and R + d are different, there will be a discontinuity in the slope 
of u n (r) at the interface and the surface integrals in Eq.s 40c and 41b will be ill-defined, depending on whether we 
use 'Pmn (r) and a mn (r) from the cell around R or R + d. Since Bloch functions do not vary greatly among different 
III-V semiconductors, we may take the surface integral to simply be a parameter depending on the atomic species at 
R and R + d. Moreover, the true u m at a heterojunction will be a smooth function depending primarily on the type 
of the dimer, with some small dependence on the nearby atoms. Taking the coefficient to depend only on the two 
atoms connected is essentially ignoring the influence of nearby atoms on the microscopic Vo ( r ) and assuming the Bloch 
function at a dimer is the same as what it would be for that dimer in a bulk binary material. The discontinuity in the 
Bloch functions at a heterojunction potentially spoils current conservation which requires ip(r) have continuous first 
derivatives. A discontinuity in the slope of u m requires a compensating discontinuity in the slope of the associated 
envelope. Such an envelope is certainly possible in the continuum, however this cannot be accomplished in a real-space 
formulation on a grid since such a discontinuity would be over a length scale smaller than the cutoff imposed by the 
grid. The Burt-Foreman formulation of envelope theory resolves this problem and ensures current conservation, but 
requires additional momentum matrix elements between wave vectors outside the first Brillouin zone[42, 43]. 


VIII. EIGHT-BAND MODEL 

We now apply the ideas developed in Secs. I-VII to the eight-band Kane model with spin-orbit coupling and 
perturbative remote band contributions. This model has been used to describe electronic states in bulk materials, 
impurities, and nanostructures. The Hamiltonian is given by[44-46] 
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where £7 C , 1%, and E a = E v — A are the zone-center energies of the conduction, valence, and split-off bands respectively, 
F is the remote band contribution to the conduction band effective mass, iPo = ^■(5 , |P a; |X), and B is the inversion 
asymmetry parameter due to remote bands, which we take to be zero. The modified Luttinger parameters are given 
by 
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Note that we use the definitions from Ref. 44 rather than Ref. 46 since they give the standard relationships between 
hole masses and Luttinger parameters, 
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For actual III-V materials the modified Luttinger parameters in Ref. 46 give effective masses that differ from the 
relationships given above by about 10%. 

To take the atomistic limit of the eight-band model we proceed as in Sect. V, replacing k’s with the appropriate 
difference operators on the crystal lattice. As in Sect. V we must include the additional atomistic momentum matrix 
element iQ a = ^ (X\p y \Y)n 1 = —^{X\p y \ Y)q 2 where Qi (ST 2 ) is the volume around the anion (cation). With this 
sign convention a positive Q a will give a long wavelength 16-band model with Q > 0 in agreement with Ref. 33. With 
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the inclusion of spin and transforming to the total angular momentum basis in which Eq. 44 is written, we obtain 
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where k± = ( k x ± ik y ). As in the case of the four-band model, Hq is the atomistic contribution on the anion site 
with —Hq on the cation site. 

Rather than using modified Voronoi cells as we did with the four-band model, here we avoid the Hermiticity problem 
by choosing a basis of inversion symmetric Bloch functions. The eight-band model with B = 0 is inversion symmetric, 
so adopting such a basis is quite natural, and in any case the choice of basis is arbitrary. Although the basis is 
symmetric, inversion symmetry is broken by Hq (c.f. Eq. la) which will cause the envelope functions to be inversion 
non-symmetric. The eight-band model depends on eight parameters (E c , E v , E a , Pq, 71 , 72 , 73 , and F). As discussed 
in Sect. V, the diagonal energies will be different on each atom, doubling their number to six. The other parameters, 
which depend on derivatives of the inversion symmetric Bloch functions, will be the same on each atom. With the 
inclusion of Q a , the number of parameters is increased from eight to 12 in the atomistic limit. 

The remote band contributions in the eight-band model introduce additional fcs, and thus additional difference 
operators. As discussed in Sect. Ill, the fact that each atom has 4 neighbors means that only fc 2 , k x , k y , and k z can 
be constructed using nearest neighbor differences while k x k y , k x k z , k y k z , fc 2 , fc 2 , and fc 2 require next nearest neighbor 
differences. Replacing the fcs with difference operators acting on plane waves as in Sect. V, the Hamiltonian becomes 
a 16 x 16 matrix of the form 


( H n H 12 \ 

V h\ 2 h 22 ) 


(49) 


where the diagonal blocks H n and H 22 include on-site and next-nearest-neighbor couplings and the H\ 2 block contains 
nearest neighbor couplings. For fc = 0 the nearest neighbor couplings from k x , k y , and k z vanish but those from fc 2 
remain, shifting the zone-center energies from their eight-band values. This can be seen in Eq. 26 where Ts are 
nonzero even for k = 0. For zincblende the onsite energies are different on the two atoms (E c 1 , E v 1 , E s \ on the anion 
and E c2 , E v2 , E s2 on the cation), breaking the inversion symmetry. With two grid sites per unit cell the number of 
bands in the original fc • p model is doubled in the atomistic limit and each zone-center state of the continuum model 
splits into two states: one having an envelope with the same sign on each atom (++) and another with opposite 
signs (H— ). The +- 1 - eigenvectors correspond to the original eight-band states (r 6c , Fy,,, r 8t ,) and the -I— states to 
correspond to the additional states found in the 16-band model (Tgi,, Tgc, r 7 C ) as shown in Fig. 2. 

To understand the band doubling, consider a continuum Hamiltonian having a diagonal term H = E c + C c k 2 where 
the subscripts indicate the parameters are for the continuum model. The parameters of the atomistic model are 
determined by fitting to the parameters of the continuum model. Because V 2 couples the two atoms within a unit 
cell the atomistic Hamiltonian is a non-diagonal 2x2 matrix even at k = 0, 
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where E a \, E a2 are the atomistic energies for atom 1 and 2 (anion and cation respectively for zincblende) and C a is the 
atomistic coefficient corresponding to C c in the continuum Hamiltonian. The zone-center eigenvalues and eigenvectors 
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FIG. 2: In the atomistic limit the number of bands is doubled due to there being two atoms per unit cell. For each 
band of the original continuum k ■ p model there will be one state with an envelope that has the same sign on both 
atoms (++) and another with an envelope having opposites signs (-1—) on the two atoms. The +-1- solutions are 
taken to be the states of the original continuum k ■ p model and the H— solutions, which are shifted in energy by 
~ 0(h 2 /m 0 a 2 att ), are taken to be excited bands. The energy of an excited band also depends on the remote band 
contribution to the effective mass (i.e. the coefficient of V 2 in the atomistic Hamiltonian). The remote contribution 
to the valence band is multiplied by y a i > 0, shifting the Tgc and Tt c states to higher energies. Because the remote 
contribution to the conduction band mass is negative the r6„ band is displaced downward. 
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In the inversion symmetric case, E a i = E a2 and the zone center energies are given by E = E a and E = E a + 
32 C a /af att . The atomistic parameter E a is then set equal to the continuum zone center energy, E c , and the atomistic 
limit gives rise to an additional band with zone-center energy 32 C a /af att higher. If C a is fixed by fitting to an effective 
mass, there are no additional parameters to fit. In the inversion non-symmetric case the zone center energies of the 
atomistic model are functions of three parameters, E a i, E a 2 , and C a , the last of which would be fixed by fitting to 
an effective mass. We require two conditions to fit E a \ and E a 2 . The simplest approach is to set Eq. 51a to E c and 
the energy of one additional band. Alternatively, one could fit to E c , and another condition such as the ratio of of 
envelope functions on the two atoms. Since empirical data for energies of excited bands is more readily available than 
Bloch functions, we will determine E a 1, and E a 2 by fitting to the two band energies. Our fitting procedure is by no 
means unique, and additional data may make some other method preferable. 

We take as our target energies those of the 16-band model, E 6v , Ej v , E Sv , E ec , E^c, E 8c (see Fig. 2). These 
have been measured for some materials [47] and calculated for others[2]. Setting the zone-center eigenvalues of the 
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atomistic Hamiltonian equal to the target energies gives the atomistic on-site energies 


Ed = E 6c + E 6v ) - 16(1 + 2 F a )h 2 /m 0 al tt - ^ (E 6c - E 6v ) 2 - 1024(1 + 2F a ) 2 /P/TO 2 af att 
E c 2 = ~{Eq c + Eq v ) — 16(1 + 2F a )H 2 /moaf att + — yj ( Eq c — Eq v ) 2 — 1024(1 + 2F a ) 2 h7 / rn^a^ att 
Esi = \(E?e + E 7v ) + 16-iaih 2 /m. 0 a 2 latt - \yj (E 7c - E 7v ) 2 - 10247^ ft 4 /mgaf att 

Es 2 = ^(-^70 + E 7v ) + 167 a ift 2 /m 0 af att + \yj (E 7c - E 7v ) 2 - 10247^ ft 4 /m;]af aM 

£,i = ^(-Esc + #8,) + /m. 0 a 2 latt - \ \J(E S , C - E Sv ) 2 - 10247^ ft 4 lm\a\ att 

E v 2 = ^ (-Bsc + ^8t>) + 167 a ift 2 / moa 2 att + I y(P 8c - P 8l ,) 2 _ 102472^4 /mlaf att (52) 

where F c i is the on-site conduction band energy on the atom at the origin (anion for zincblende), E c 2 is the on-site 
energy for the atom at (cq aM /4, a/ ott /4, a; att /4) (cation for zincblende), and likewise for the valence band (subscripts 
id and v2) and the spin-orbit band (subscripts si and s2). The a subscripts on F„ and 7 a i indicate they are the 
atomistic versions of the continuum k ■ p parameters F and 71 . The corresponding envelope functions are obtained 
from Eq. 51b, giving 

die, = ([E c2 - Ed)m 0 a 2 att + 1024(1 + 2F a ) 2 ft 4 + ( E c2 - F cl ) 2 m 2 a 4 at ^ /32(1 + 2 F a )H 2 

A 6c = ({E c 2 ^ Ed)m 0 al tt - 1024(1 + 2F a ) 2 /H + (F c2 - F cl ) 2 m 2 a 4 att ) /32(1 + 2 F a )H 2 

A 7v = - (^(E s 2 - E sl )m 0 af att + yj \02^ 2 al tA + ( E s2 - E sl ) 2 mlaf at jj /32 7a i h 2 

A 7c = - ^( E s2 - E sl )m 0 af att - yJ\02^ 2 al h A + ( E s2 - E sl ) 2 mlaf at jj /32^ al h 2 

A 8v = - (^{E v2 - E vl )m 0 af att + yj 10247 ^h 4 + ( E v2 - E vl ) 2 mlaf at jj /32j al h 2 

A Sc = - ([E v2 - E vl )m 0 a 2 att - yj 10247^/i 4 + (E v2 - E vl ) 2 mlaf at jj /32'y al H 2 . (53) 

The states with energies Eq c , E 7v , E 8v have envelopes of the form +- 1 - while those with energies Eq v , F 8c , E 7c have 
H— envelopes. Note that the envelopes break inversion symmetry and the momentum matrix elements taken over a 
unit cell must include the sub-unit cell structure of the envelopes. 

Having determined the on-site energies, we must now determine the coefficients of the difference operators. In k ■ p 
theory effective masses only require knowing E{k) to second order in k , and are therefore computed using second 
order perturbation theory in A: [34, 48, 49]. The g-factor is also computed this way [49]. For the 16-band continuum 
k ■ p model the results are 


^ = (i + 2F) + % 
m* 6 

9* , . 9r E Po 

3 


( 2 1 1 ) 

\ + Epi \ 

( 2 1 1 ^ 

V Eq c Eg v E 6c Ej v J 

l+ 3 1 

^ Eq c — e 8c Eq c — e 7c J 


9o 

7i = 7i 


72 = 72 


73 = 73 


9o 


( 1 1 ) 

1 E Pl 

f 1 1 ^ 

V Eqc — e 8v E 6c — e 7v ) 

] 3 1 

^ Eq c — e 8c e 8c — e 7c J 


Ep 0 


Ep 2 


Eq 


E, 


Q 


3(E$ V — Eq c ) 3 (E$ v ~ Eq v ) 3 {E 8v — E 7c ) 3 (E 8v — E 8c ) 

Ep 0 _ Ep 2 Eq 

6(Esv — Eq c ) 6(F 8 „ — Eq v ) 6(F 8 „ — E 7c ) 

Epp _ Ep 2 _ Eq 

6(F 8 „ — Eq c ) 6(F 8 „ — Eq v ) 6 (E 8v — E 7c ) 


(54a) 

(54b) 

(54c) 

(54d) 

(54e) 


where go is the bare electron g factor, g r is a possible remote band contribution, Ep 0 = 2|(S' c |P r |X u )| 2 /TOo, Ep 1 = 
2|(5 , c |P x |X c )| 2 /m 0 , Ep 2 = 2|(5' t ,|P a; |X u )|2/m 0 , and Eq = 2|(X c |P y |Z„)| 2 /mo. The quantities on the left hand sides 
of Eq.s 54a-54e are the target parameters taken from experiment or possibly ab initio calculations while P, Ep 0 , 
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Ep 1 , Eq , 71 , 72 , and 73 are the model parameters empirically chosen to reproduce the target values. The model 
parameters are easily determined since the equations are linear. Since E 8v — Eq v is large, we may set Ep 2 = 0, giving 
five equations in five unknowns. 

The effective masses and g-factors of the atomistic model may also be computed using perturbation theory for small 
k. The resulting expressions contain effective matrix elements that depend on the variation of the envelope over the 
unit cell and the bands which they connect. For example, the momentum matrix element between the Tg c and Tg,, 
states is 

P6c8v = ~i — (£r 6c |-Px|^r 8 ,,} 
mo 

1 f A 6c \ T f 0 iP a \ 1 f As v \ 

1 ) 0 J yiT^V 1 ) 

= Pa(A &c + A 8v )/yj (1 + A\ c ){\ + A\ v ) (55) 


with the other matrix elements {Pq C 8v, Pq C 7 v, PecSc, Pqc 7 c, P&v8v, Pdv 7 v , Pqv8v, and P 6v 7v ) defined similarly. The Q 
matrix elements take the form 


n ■ ^ /v id 17 \ (^8 c A 8v )Q a 

Qscsv = -i—\ A r 8 c E y Zr 8 J = , ■ 

Doing second order perturbation theory directly on the atomistic Hamiltonian we obtain 
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^4-8 vQ a^Latt 


6(F'8x — E 6v ) 

Ep 8v 6c 
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Ep 8u Qv Eq g v 7 C 
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6(1 + ^4^) Q{Esv — E 6c ) 3 (E Sv — E ev ) 6 (E 8v — E 7c ) 6 (E Sv — E 7v ) 


(56) 

(57a) 

(57b) 

(57c) 

(57d) 

(57e) 


where the band-specific Kane energies are Ep mn = 2\{Sm\Px\X n ) | 2 /too = 2moP% ln /h 2 and EQ mn = 
2\{X m \p y \Z n )\ 2 /m.Q = 2moQ^ nn /h 2 . The perturbative expression are the same as in continuum k ■ p except that 
matrix elements are replaced with the effective matrix elements including variation of the envelope within the unit 
cell. Unlike the continuum case, the dependence on F a and 7 a i is nonlinear. 


IX. PARAMETER FITTING 

To determine the atomistic parameters we adopt a set of target material parameters to which we fit the atomistic 
parameters (see Table I). These values are known with varying degrees of certainty, with some taken from high 
precision measurements while others are obtained theoretically. The basic eight-band parameters are taken from Ref. 
50. The zone-center energies of the higher lying bands (Tg,,, T 7 C , Tg c ) are not as well known and we have taken their 
values from the k ■ p calculation of Ref. 6 , with the exception of GaAs for which we used the experimental values 
from Ref. 47. The values of the conduction band effective g factor, g *, and the Dresselhaus spin splitting, 7 C , were 
also taken from Ref. 6 . Some modifications have been made for InAs. The spin orbit coupling has been increased to 
A = 0.45 eV in order to be able to obtain g* = —14.9 without having 7 C > 100 eV A 3 . Alternatively, the A could be 
left unchanged and g* fit using a remote band contribution g r . The value of 7 ^ for InAs has been reduced from 8.2 
to 7.5 in order to avoid bands that cross the gap at large fc, which is much less of a liberty than it may seem since 
the Luttinger parameters for InAs are poorly known[50]. 

The atomistic on-site potentials (£ c i, E c 2 , E v 1 , etc. in Table II) are determined from Eq. 52 using the lattice 
constants and zone-center energies from Table I. To fit the effective masses we must determine the values of F a , P a , 
Q a , 7 a 1 j 7 a 2 ) 7 a 3 for which Eq. 57a-57e match the empirical target values. This is more difficult than the fitting 
procedure for a continuum k-p model because the effective momentum matrix elements and remote band contributions 
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depend on F a and 7 a i. We do a nonlinear fit on F a and 7 a i. For particular values of F a and 7 a i, P a is determined by 
Eq. 57a, Q a is determined by Eq. 57c, and y a 2 and 7 a 3 are determined by Eq.s 57d and Eq. 57e. F a is then adjusted 
to make the resulting g* match the target value. This results in a curve in the F a - 7 al plane from which we pick the 
point at which the Dresselhaus spin splitting fits the target as well. We determine y c by numerically computing the 
spin-splitting E(k) in the 110 direction. The range of F a and 7 a i which must be numerically searched is reduced 
by the condition that the amplitudes A § v ,..., As c must be real and the solutions corresponding to E(T6 C ), E(T 7v ), 
E(T Sv ) must have signature +-(-. These conditions restrict the values to —af att (E 6c — E 6v )/32 < (2 F a + 1) < 0 and 
~ a latt(E'8c ~ Eg v )/32 < 7 a i < 0. 

The band structures resulting from our numerical fits are shown in Fig. 3. Since the atomistic model is derived 
from a continuum k ■ p model that is perturbative in k, our results are accurate for small k and all materials appear 
to have a direct gap. Perturbative k ■ p models eventually break down at large k, which can result in spurious 
solutions that cross the gap. When working in a plane wave basis these spurious solutions may be avoided by simply 
restricting the values of k, however this cannot be done in a real-space formulation. Spurious gap-crossing states can 
be eliminated by modifying the basis [51, 52], choosing different material parameters [5 ], or altering the differencing 
scheme to include higher powers of k that push the spurious solutions out of the gap[54]. The threat of gap-crossing 
bands is greater in the atomistic limit due to the larger ( computational ) Brillouin zone associated with the smaller 
computational grid, providing more space for the bands to turn over and cross the gap. We show energies throughout 
the entire Brillouin zone to demonstrate that our parameterization does not produce spurious gap-crossing states, 
and the model is suitable for use in a real-space formulation. Only InAs required modifications to the parameters to 
suppress spurious solutions, as described above. 


parameter 

A1P 

GaP 

InP 

AlAs 

GaAs 

InAs 

AlSb 

GaSb 

InSb 

ai a tt (A) “ 

5.4584 

5.4417 

5.8613 

5.6524 

5.6416 

6.0501 

6.1277 

6.0817 

6.4690 

E 6v (eV) 6 

-11.21 

-12.14 

-11.04 

-11.73 

-12.9 C 

-11.53 

-10.62 

-11.47 

-10.54 

E 7v (eV) “ 

-0.07 

-0.08 

-0.108 

-0.28 

-0.341 

-0.45 d 

-0.676 

-0.76 

-0.81 

E 8v (eV) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

E 6c (eV) “ 

3.63 

2.886 

1.423 

3.099 

1.519 

0.417 

2.386 

0.812 

0.235 

E 7c (eV) e 

4.78 

4.38 

4.78 

4.55 

4.488 c 

4.858 d 

3.53 

3.11 

3.18 

E 8 c (eV) e 

4.82 

4.47 

4.97 

4.70 

4.659 c 

4.79 d 

3.77 

3.44 

3.64 

m*/m 0 (1) a 

0.22 

0.13 

0.0795 

0.15 

0.067 

0.026 

0.14 

0.039 

0.0135 

g*/go (i) e 

1.92 

1.9 

1.26 

1.52 

-0.44 

-14.9 

0.84 

-9.2 

-51.6 

ji (i) “ 

3.35 

4.05 

5.08 

3.76 

6.98 

20 

5.18 

13.4 

34.8 

72 (1) “ 

0.714 

0.49 

1.6 

0.82 

2.06 

7.5 / 

1.19 

4.7 

15.5 

73 L (1) “ 

1.23 

1.25 

2.1 

1.42 

2.93 

9.2 

1.97 

6.0 

16.5 

7c (eVA 3 ) “ 

2.1 

-2.4 

-8.4 

11.4 

25.0 

40.5 

40.9 

185.0 

226.0 


“Ref. 50, except where noted for InAs. 

^Ref. 55, except for InAs. 

“Ref. 47. 

^Modified to fit g* and 7 C . 

“Ref. 6 except where noted for GaAs and InAs. 

^Modified to avoid gap-crossing bands at large k. 

TABLE I: Target material parameters. These are the physical material parameters which the atomistic model 
parameters were adjusted to fit. The values of g*/go and y c from Ref. 6 are experimental values when available, and 

theoretical values when no measurements were available. 


X. CONCLUSION 

We have demonstrated how to construct an atomistic k ■ p theory with finite differences on a grid matched to 
the crystal lattice. Taking the atomistic limit of k ■ p theory in a straight-forward way results in a non-Hermitian 
Hamiltonian, which is seen to be related to the well known fact that multiplying a difference operator by a spatially 
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parameter 

A1P 

GaP 

InP 

AlAs 

GaAs 

InAs 

AlSb 

GaSb 

InSb 

E cl (eV) 

1.560633 

0.516670 

-1.501874 

1.387058 

0.675926 

-0.354316 

-1.094007 

-0.568950 

-11.520260 

E c2 (eV) 

5.236823 

4.669935 

3.360307 

4.483957 

2.273318 

1.099526 

4.568593 

1.934999 

-0.968967 

E v i (eV) 

-0.505855 

-0.593772 

-0.604112 

-0.302607 

-0.751431 

-0.308368 

-0.328926 

-0.434964 

-0.310302 

E v 2 (eV) 

0.645350 

0.824859 

0.809339 

0.347832 

1.162880 

0.353026 

0.400082 

0.591758 

3.264512 

E s i (eV) 

-0.620739 

-0.597985 

-0.591704 

-0.836975 

-1.198921 

-0.788711 

-1.423545 

-1.490988 

-1.213870 

E s2 (eV) 

0.650234 

0.689072 

0.498932 

0.452200 

1.098370 

-0.056631 

0.578701 

0.557782 

2.808081 

Fa (1) 

-1.378388 

-1.376856 

-1.308404 

-1.450082 

-1.435248 

-1.390038 

-1.401502 

-1.411955 

-0.312572 

7 al (1) 

-0.571909 

-0.514787 

-0.671323 

-0.609904 

-0.554424 

-0.731854 

-0.569585 

-0.498025 

-0.117695 

7^2 (1) 

0.118013 

-0.733646 

-0.582392 

-0.292962 

-0.607743 

0.572230 

0.201367 

-0.087151 

-1.036618 

7a3 (1) 

0.139473 

-0.334737 

-0.221699 

-0.010489 

-0.113528 

0.633312 

0.260787 

0.246206 

-0.582445 

P a (eVA) 

9.325446 

10.188130 

8.866846 

10.192099 

10.396101 

8.791500 

9.428844 

10.123984 

9.703173 

Qa (eVA) 

-7.062907 

-5.962138 

-5.085142 

-6.279963 

-6.139187 

-10.995692 

-7.559866 

-7.470334 

-4.704146 


TABLE II: Empirical atomistic k.p theory parameters fit to the target parameters in Table I 


varying coefficient leads to non-Hermiticity. A more careful treatment shows the problem may be remedied by using 
the finite volume method, starting with an inversion symmetric Bloch basis, or by deforming the computational cells to 
generalized Voronoi cells. The use of symmetric Bloch functions does not limit one to the symmetric approximation 
since the atomistic envelope functions themselves vary within the unit cell even at k = 0. The use of inversion 
symmetric Bloch functions and generalized Voronoi cells solve the Hermiticity problem, but are not applicable to 
heterojunctions. As a result these approaches can be used on systems such as bulk materials, bulk materials with 
impurities or applied potentials, or nanocrystals with a vacuum barrier. The finite volume method can be used in the 
presence of heteroj unctions. 

The atomistic limit of a simple four-band k-p model exactly reproduces the four-band tight-binding model, provided 
we include spherically symmetric remote band contributions for both the conduction and valence band, and the 
atomistic momentum matrix elements are different on different atoms (for zincblende). In order to have different 
momentum matrix elements that make the model exactly match the tight-binding model requires the use of generalized 
Voronoi cells to to symmetrize the momentum matrix elements without making them all equal. The atomistic limit 
of the widely used eight-band model results reproduces effective masses, g-factors, and Dresselhaus spin splittings of 
III-V materials. The fits are exact for most materials, with the exception of InAs for which it was necessary to increase 
the spin-orbit coupling. This may be due to an insufficient number of bands in the model, or due to uncertainties in 
the experimental values. 

The particular implementation presented in Sect. VIII is by no means unique, and different atomistic models are 
possible depending on the choice of Bloch basis (inversion symmetric or not), the differencing scheme, and whether 
or not remote band contributions are included. In addition, different fitting procedures may be used. For example, 
the higher lying band energies could be left as free parameters adjusted to fit the band structure to other criteria 
such as charge asymmetry. We have chosen to exactly fit all zone center energies, the zone center effective masses 
for the bottom of the conduction and top of the valence bands, as well as conduction g-factors and Dresselhaus spin 
splittings, since these quantities are the most important for electronic states of impurities and nanostructures. 

An interesting property of these models is that the envelope functions have momenta outside the first Brillouin 
zone, a feature shared with the Burt-Foreman[42, 43] approach to dealing with heterojunctions. Since our model is 
constructed in real space there is not a clearly defined separation between the wave function components that are 
associated with Bloch functions and those that are not, while the Burt-Foreman approach has a clear distinction 
between Bloch and envelope functions in fc-space. 

As seen from the fit to III-V materials, the atomistic envelope theory can reproduce the effective masses of the bands 
near the gap. This is in contrast to tight-binding models which can give incorrect effective masses[56]. The four-band 
model with only spherically symmetric remote band contributions illustrates how a nearest neighbor tight-binding 
model fails to reproduce the correct cubic band warping of the valence band. In contrast, the atomistic Kane model 
gives the correct effective masses because it contains next nearest neighbor couplings via the Luttinger parameters. 

There are many potential applications of this method to the electronic properties of impurity states, alloys, and 
polytypes[57]. For sufficiently small nanoparticles we expect atomistic k ■ p theory to improve the description of the 
electronic structure, compared with continuum k ■ p-theory. In particular nanoparticles with an irregular surface 
necessitate an atomistic description. It would be interesting to test how well our new method can describe structural 
defects such as dislocations, twin planes and stacking defects. Such systems cannot be easily treated in continuum 
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FIG. 3: Band diagrams using the parameters in Table II. Because the atomistic model is derived from a model that 
is perturbative in k, it is inaccurate for large k and spurious solutions can cross the gap. In a real-space formulation 
such spurious solutions will contribute to, and thus spoil, a numerical solution. It is therefore important to examine 
the entire Brillouin zone associated with the computational grid (or the crystal lattice in the atomistic limit). No 
gap-crossing states are seen for the parameters given in Table II. InAs required an adjustment of the material 

parameters to avoid gap-crossing states (see Sect. IX ). 


k.p-theory. Atomistic k ■ p theory will also allow strain effects to be directly modelled in terms of atomic positions, a 
task which is difficult in both tight-binding and pseudopotential methods. Finally, atomistic k-p theory has the unique 
feature that it allows the combination of atomistic and continuum models in the same system to facilitate multiscale 
modeling since the grid can be highly non-uniform. One could use a rectilinear grid in ’’large” regions described by 
a continuum model and an atomistic grid in the regions requiring atomistic precision. The differencing operators in 
the regions where the rectilinear and atomistic grids meet would be peculiar to the details of the grid used, but would 
be well defined. Multiscale modeling will dramatically reduce the computing time of atomistic k.p-theory compared 
with other atomistic models, while keeping atomistic accuracy where it is necessary. 

Acknowledgements M.-E. P. acknowledges the support of the Swedish Research Council (VR). 
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XI. APPENDIX: BLOCH BASIS STATES 


The Bloch state basis for the eight-band model is 


Ml = M r _« /2 = \S 4.) 

«2 = W+i/2 = | S t) 



u 6 =u _\ /2 = -j=\ (X-iY) t)+i\J-\Zl) 
u 7 = u r _\ /2 = -J|(X - iY) t) + ^\Z i) 
u s = < 7 1/2 = ^ | (X + iY) i) - -±= \Z f) 


\Zt) 


where the ordering of states is the same as for the Hamiltonian in Eq. 44. 
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